###Full-Power Cabinets and Caretaker Administrations in Parliamentary Democracies, 1945-2024###
###Francesco Bromo and José A. Cheibub###
###15 July 2025###

rm(list = ls())

##Load packages
library(haven)
library(plm)
library(dplyr)
library(psych)
library(mgcv)
library(ggplot2)
library(ggpubr)

##Import data
formation_data <- read_dta("formation_data.dta") 

#Table A3: Summary Statistics 
describe(formation_data$formation_days) #Caretaker Duration (Days)
describe(formation_data$investiture) #Investiture
describe(formation_data$post_elec) #Post-Election
describe(formation_data$enp) #ENPP
describe(formation_data$sd_rile) #Ideological Polarization 

##Generalized Additive Models
data <- pdata.frame(formation_data, index=c("country_number", "event_date")) 

data$id_year <- data %>% group_indices(event_year) 

data$country_number[data$country=="Australia"] <- 31
data$country_number[data$country=="Slovenia"] <- 1

#Table A4, column 1 (Whole Sample)
gam1 <- gam(formation_days ~ investiture + post_elec + enp +  sd_rile + enp*post_elec + sd_rile*post_elec + s(id_year, bs="cr") + factor(country_number), data = data)
summary(gam1)

#Table A4, column 2 (Established Democracies)
data_estd <- data[which(data$europe != 2),]

gam2 <- gam(formation_days ~ investiture + post_elec + enp +  sd_rile + enp*post_elec + sd_rile*post_elec + s(id_year, bs="cr") + factor(country_number), data = data_estd)
summary(gam2)

#Table A4, column 3 (Eastern Europe)
data_ee<-data[which(data$europe==2),]

gam3 <- gam(formation_days ~ investiture + post_elec + enp +  sd_rile + enp*post_elec + sd_rile*post_elec + s(id_year, bs="cr") + factor(country_number), data = data_ee)
summary(gam3)

#Table A5: Country Dummies (Table A4)
coefs <- as.data.frame(summary(gam1)$p.table)
country_coefs <- coefs[grep("^factor\\(country_number\\)", rownames(coefs)), ]
country_numbers <- as.numeric(gsub("factor\\(country_number\\)", "", rownames(country_coefs)))

country_map <- unique(data[, c("country_number", "country")])

country_dummies <- data.frame(
  country_number = country_numbers,
  estimate = country_coefs[, "Estimate"],
  std_error = country_coefs[, "Std. Error"]
)

merge(country_dummies, country_map, by = "country_number") %>%
  select(country, estimate, std_error)

#Figure A1, panel A
p_obj <- plot(gam1, residuals = TRUE)
p_obj <- p_obj[[1]]
sm_df <- as.data.frame(p_obj[c("x", "se", "fit")])

plot_whole <- ggplot(sm_df, aes(x = x, y = fit)) +
  geom_ribbon(aes(ymin = fit - se, ymax = fit + se, y = NULL), alpha = 0.125) +
  geom_line(size = 1) +
  labs(x = "Time", y = "s(Time)") +
  theme_classic2() +
  theme(axis.title = element_text(size = 15, color = "black"), axis.text = element_text(size = 15, color = "black")) +
  ggtitle("A. Whole Sample") +
  theme(plot.title = element_text(size = 15, face = "bold", hjust = 0.5)) +
  coord_cartesian(ylim = c(-25, 30)) + scale_x_continuous(breaks = seq(0, 75, 15), labels = c("0" = "1945", "15" = "1960", "30" = "1975", "45" = "1990", "60" = "2005", "75" = "2020"))
plot_whole

#Figure A1, panel B
p_obj_estd <- plot(gam2, residuals = TRUE)
p_obj_estd <- p_obj_estd[[1]]
sm_df_estd <- as.data.frame(p_obj_estd[c("x", "se", "fit")])

plot_estd <- ggplot(sm_df_estd, aes(x = x, y = fit)) +
  geom_ribbon(aes(ymin = fit - se, ymax = fit + se, y = NULL), alpha = 0.125) +
  geom_line(size = 1) +
  labs(x = "Time", y = "s(Time)") +
  theme_classic2() +
  theme(axis.title = element_text(size = 15, color = "black"), axis.text = element_text(size = 15, color = "black")) +
  ggtitle("B. Established Democracies") +
  theme(plot.title = element_text(size = 15, face = "bold", hjust = 0.5)) +
  coord_cartesian(ylim = c(-25, 30)) + scale_x_continuous(breaks = seq(0, 75, 15), labels = c("0" = "1945", "15" = "1960", "30" = "1975", "45" = "1990", "60" = "2005", "75" = "2020"))
plot_estd

#Figure A1, panel C
p_obj_ee <- plot(gam3, residuals = TRUE)
p_obj_ee <- p_obj_ee[[1]]
sm_df_ee <- as.data.frame(p_obj_ee[c("x", "se", "fit")])

plot_ee <- ggplot(sm_df_ee, aes(x = x, y = fit)) +
  geom_ribbon(aes(ymin = fit - se, ymax = fit + se, y = NULL), alpha = 0.125) +
  geom_line(size = 1) +
  labs(x = "Time", y = "s(Time)") +
  theme_classic2() +
  theme(axis.title = element_text(size = 15, color = "black"), axis.text = element_text(size = 15, color = "black")) +
  ggtitle("C. Eastern Europe") +
  theme(plot.title = element_text(size = 15, face = "bold", hjust = 0.5)) +
  coord_cartesian(ylim = c(-25, 30), xlim = c(0, 75)) + scale_x_continuous(breaks = seq(0, 75, 15), labels = c("0" = "1945", "15" = "1960", "30" = "1975", "45" = "1990", "60" = "2005", "75" = "2020"))
plot_ee

#Figure A1: Effect of Time on the Duration of Caretaker Administrations, 1945-2020
plot_all <- ggarrange(plot_whole, plot_estd, plot_ee, ncol = 2, nrow = 2)
plot_all

#getwd()
#ggsave("plot_all.jpeg", plot_all, width=7.5, height=5, dpi=300)